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A density functional theory for many-body lattice models is considered in which the single-particle 
density matrix 7ij is the basic variable. Eigenvalue equations are derived for solving Levy's con- 
strained search of the interaction energy functional W^[7ij]. VK[7ij] is expressed as the sum of 
Hartree-Fock energy Eupl'Jij] and the correlation energy Ec['yij]- Exact results are obtained for 
£0(712) of the Hubbard model on various periodic lattices, where yij — 712 for all nearest neigh- 
bors i and j. The functional dependence of £0(712) is analyzed by varying the number of sites 
Na, band filling A'^e and lattice structure. The infinite one-dimensional chain and one-, two-, or 
three-dimensional finite clusters with periodic boundary conditions are considered. The properties 
of £0(712) are discussed in the limits of weak (712 ~ 7^2) and strong (712 ~ ■y^) electronic correla- 
tions, and in the crossover region (712 < 712 < 712)- Using an appropriate scaling we observe that 
£0(512) ~ Ec/Ehf has a pseudo-universal behavior as a function of (712 = (712 — 7i2)/(7i2 ~ Til)- 
The fact that £0(312) depends weakly on Na, N^ and lattice structure suggests that the correla- 
tion energy of extended systems could be obtained quite accurately from finite cluster calculations. 
Finally, the behavior of £0(712) for repulsive (L'' > 0) and attractive {U < 0) interactions are 
contrasted. 
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I. INTRODUCTION 

Density functional theory (DFT) has been the subject 
of remarkable developments since its original formulation 
by Hohenberg and Kohn (HK) ||l|,0|. After formal im- 
provements, extensions, and an uncountable number of 
applications to a wide variety of physical problems, this 
theoretical approach has become the most efficient, al- 
beit not infallible, method of determining the electronic 
properties of matter from first principles 0,0] . The most 
important innovation of DFT, which is actually at the 
origin of its breakthrough, is to replace the wave func- 
tion by the electronic density p{r) as the fundamental 
variable of the many-body problem. In practice, den- 
sity functional (DP) calculations are largely based on the 
Kohn-Sham (KS) scheme that reduces the many-body iV- 
particle problem to the solution of a set of self-consistent 
single-particle equations ||]. Although this transforma- 
tion is formally exact, the implementations always re- 
quire approximations, since the KS equations involve 
functional derivatives of the unknown interaction energy 
14^[p(r)], usually expressed in terms of the exchange and 
correlation (XC) energy Excipi"!^)]- Therefore, under- 
standing the functional dependence of i?xc[p('^)] a-nd im- 
proving its approximations are central to the develop- 
ment of DP methods. The currently most widespread 
Ansdtze for Exc [p(^)] — the local density approximation 
(LDA) with spin polarized [|| and gradient corrected 
extensions p] — were originally derived from exact re- 
sults for the homogeneous electron gas. It is one of the 



purposes of this paper to investigate the properties of the 
interaction-energy functional from an intrinsically inho- 
mogeneous point of view, namely, by considering exactly 
solvable many-body lattice models. 

Despite the remarkable success of the local spin den- 
sity approximation, present DPT fails systematically in 
accounting for phenomena where strong electron correla- 
tions play a central role, for example, in heavy-fermion 
materials or high-Tc superconductors. These systems 
are usually described by simplifying the low-energy elec- 
tron dynamics using parameterized lattice models such 
as Pariser-Parr-Pople, |7| Hubbard, g or Anderson g] 
models and related Hamiltonians [Q . Being in principle 
an exact theory, the limitations of the DP approach have 
to be ascribed to the approximations used for exchange 
and correlation and not to the underlying HKS formal- 
ism. It would be therefore very interesting to extend the 
range of applicability of DPT to strongly correlated sys- 
tems and to characterize the properties E'xc in the limit 
of strong correlations. Studies of the XC functional on 
simple models should provide useful insights for future 
extensions to realistic Hamiltonians. Moreover, taking 
into account the demonstrated power of the DP approach 
in ah initio calculations, one may also expect that a DPT 
with an appropriate Exc could become an efficient tool 
for studying many-body models, a subject of theoretical 
interest on its own. 

Several properties of DPT on lattice models have been 
already studied in previous works 0-|l3]. Gunnarsson 
and Schonhammer were, to our knowledge, the first to 



propose a DF approach on a semiconductor model in or- 
der to study the band-gap problem [Q. In this case the 
local site occupancies were treated as the basic variables. 
Some years later Schindlmayr and Godby [|l2| provided 
a different formulation of DFT on a lattice by consider- 
ing as basic variables both diagonal elements ju and off- 
diagonal elements "fg of the single-particle density ma- 
trix (see also |14-16|). Schonhammer et al. then derived 
a more general framework that unifies the two previous 
approaches ||l^. Using Levy's constrained search method 
pTf they showed that different basic variables and differ- 
ent W functionals can be considered depending on the 
type of model or perturbation under study. Site occu- 
pations alone may be used as basic variables, if only the 
orbital energies are varied (i.e., if all hopping integrals 
tij are kept constant for i ^ j). However, off-diagonal 
elements of the single-particle density matrix must be 
included explicitly if the functional W is intended to be 
applied to more general situations involving different val- 
ues of tij^ for example, the Hubbard model on various 
lattice structures or for different interaction regimes, i.e., 
different U/t. 

In this paper we investigate the properties of Levy's 
interaction-energy functional 14^ as a function of 7ij by 
solving the constrained search minimization problem ex- 
actly. In Sec. P the basic formalism of density-matrix 
functional theory (DMFT) on lattice models is recalled 
and the equations for determining Vl^[7ij] are derived. 
Sec. [H presents and discusses exact results for the cor- 
relation energy Eq of the Hubbard model, which is given 
by the difference between W and the Hartree-Fock en- 
ergy -Ehf- These are obtained, either numerically for 
finite clusters with different lattice structures, or from 
the Bethc-Ansatz solution for the one-dimensional chain. 
Finally, Sec. |^ summarizes our conclusions and points 
out some relevant extensions. 



II. THEORY 



H 



In Sec. [I A the main results of Levy's formulation of 
DMFT are presented in a form that is appropriate for 
the study of model Hamiltonians such as the Hubbard 
model. Here, the hopping integrals tij between sites 
(or orbitals) i and j play the role given in conventional 
DFT to the external potential Vext{'r)- Consequently, 
the single-particle density matrix 7^- repla ces t he den- 
sity p{r) as basic variable |12-nq]. In Sec. II B, we de- 



rive equations that allow to determine Levy's interaction- 
energy functional W^[7y] in terms of the ground-state en- 
ergy of a many-body Hamiltonian with effective hopping 
integrals Ay that depend implicitly on 7jj . 



A. DMFT of lattice models 



We consider the many-body Hamiltonian 



ija 






(2.1) 



where cj^ (cia) is the usual creation (annihilation) opera- 
tor for an electron with spin a at site (or orbital) i. H can 
be regarded as the second quantization of Schrodinger's 
equation on a basis |lq ]. However, in the present pa- 
per, the hopping integrals tij and the interaction matrix 
elements Vijki are taken as parameters to be varied in- 
dependently. The matrix tij defines the lattice (e.g., one 
dimensional chains, square or triangular two-dimensional 
lattices) and the range of single-particle interactions (e.g., 
up to first or second neighbors). From the ab initio per- 
spective tij is given by the external potential and by the 
choice of the basis |18]. Vijki defines the type of many- 
body interactions which may be repulsive (Coulomb like) 
or attractive (in order to simulate electronic pairing) and 
which are usually approximated as short ranged (e.g., 
intra-atomic). Eq. (2.1) is mainly used in this section to 



derive general results which can then be applied to var- 
ious specific models by simplifying the interactions. A 
particularly relevant example, to be considered in some 
detail in Sec. |l|, is the single-band Hubbard model with 
nearest neighbor (NN) hoppings pj, which can be ob- 



tained from Eq. (2.1) by setting tij = — t for i and j 



NN's, tij = otherwise, and Vijki = US,j5kiS.,k Il7|,|l0[. 

In order to apply DMFT to model Hamiltonians of 
the form ( ^.1[ ) we follow Levy's constrained search pro- 
cedure |1^ as proposed by Schindlmayr and Godby [n2| . 
The ground-state energy is determined by minimizing the 
functional 



E[j,j]^EKhj] + W[j,j] 



(2.2) 



with respect to the single-particle density matrix jij. 
Eljij] is physically defined for all density matrices that 
can be written as 



7^j=^7y. =^(*|cIcj<t|*) 



(2.3) 



for all i and j, where |$) is an A'^-particle state. In other 
words, jij must derive from a physical state. It is then 
said to be p ure- state A^-representable |1^,^. The first 
term in Eq. (|2.2|) is given by 



E 



K 



= E^ 



ij lij 



(2.4) 



It includes all single-particle contributions and is usually 
regarded as the kinetic energy associated wit h th e elec- 
tronic motion in the lattice. Notice that Eq. (2.4) yields 
the exact kinetic energy for a given 7.^. There are no 
corrections on Ex to be included in other parts of the 
functional as in the KS approach. The second term in 
Eq. (2.2) is the interaction-energy functional given by ll^] 



W[^,,] = ' 



iT^Vr 



nrakl 



(*[7.: 



-t -t 
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I* [7..]) 



(2.5) 



The minimization in 
over all A^-particles 



Eq. (2.5) implies a search 
states r^[7ij]) that satisfy 



{^[li]]\ Y^a^ic^^J" l*[7ij]) = lij for all i and j. There- 
fore, Vl^[7ij] represents the minimum value of the inter- 
action energy compatible with a given density matrix 
7ij . W is usually expressed in terms of the Hartree-Fock 
energy 

Euvili]] = 9 X! ^^J'=' (7yrT7fc;(T' - 5aa'lil<ylkja) (2.6) 



and the correlation energy i?c[7ii] as 
VF[7»,]=i?HF[7»,]+i?cbj] 



(2.7) 



W and Eq are universal functionals of 7^ in the sense 
that they are independent of i^ , i.e., of the system un- 
der study. They depend on the considered interactions or 
model, as defined by Vijki , on the number of electrons iVe, 
and on the structure of the many-body Hilbert space, as 
given by N^. and the number of orbitals or sites Na- No- 
tice that Eq in Eq. ( p.7[ ) does not include any exchange 
contributions. Given 7^ (7^0- = ^ij/'2- in nonmagnetic 
cases) there is no need to approximate the exchange term , 
which is taken into account exactly by E-n-p [Eq. (2.6)]. 
Nevertheless, if useful in practice, it is of course possible 
to split W in the Hartree energy £'h and the exchange 
and correlation energy E'xc is a similar way as in the KS 
approach. 

The variational principle results from the following two 
relations lITTtl: 



Egs <^UJl^:j+W[-^^:j] 



for all pure-state 7V-representable 7^ [^ , and 



Er, 



i] 



unt; + w[t. 



(2.8) 



(2.9) 



^ gs\H\^ gs) refers to the ground-state en- 
\^ gs) to the ground-state 



io-'-JO- 



where Egg — 

ergyand7f/ = (^g^lX;^ 
single-particle density matrix. 

As already pointed out in previous works [^,^, W 
and Eq depend in general on both diagonal elements 
7ii and off-diagonal elements 7^ of the density-matrix, 
since the hopping integrals tij are non local in the sites. 
The situation is similar to the DF approach proposed by 
Gilbert for the study of non-local potentials Vext{r^r') 
as those appearing in the theory of pseudo-potentials 
|0-|l6| . A formulation of DFT on a lattice only in terms 
of 7ii would be possible if one would restrict oneself to 



a family of models with constant tij for i ^ j ■ However, 
in this case the functional Wjji,] would depend on the 
actual value of ty for i ^ j [13[ . 

The functional M^[7y], valid for all lattice struc- 
tures and for all types of hybridizations, can be simpli- 
fied at the expense of universality if the hopping inte- 
grals are short ranged. For example, if only NN hop- 
pings are considered, the kinetic energy Ek is inde- 
pendent of the density-matrix elements between sites 
that are not NN's. Therefore, the constrained search 



in Eq. (2.5) may restricted to the |^l'[7y]) that satisfy 

(*[7u]| E^Cci^T |*[7ij]> = lij only for i = j and for 
NN ij. In this way the number of variables in Vl^[7ij] is 
reduced significantly rendering the interpretation of the 
functional dependence far simpler. While this is a great 
practical advantage, it also implies that W and Ec lose 
their universal character since the dependence on the NN 
7ij is now different for different lattices. In Sec. Ill re- 
sults for one-, two-, and three-dimensional lattices with 
NN hoppings are compared in order to quantify this ef- 
fect. 



For the applications in Sec. HI we shall consider the 



single-band Hubbard model with NN hoppings, which in 
the usual notation is given by B 



{i,J)o- i 



(2.10) 



In this case the interaction energy functional reads 



W[i^,] = 



VY.^^^^'^.Wni^ni^ |*[7.j]) 



, (2.11) 



where the minimization is performed with respect to all 
TV-particle |*[7„-]) satisfying (*[7zj]| E, cLcJ<.|*[7^J])^ = 
7y- for i and j NN's. If the interactions are repulsive 
([/ > 0) M^[7y] represents the minimum average num- 
ber of double occupations which can be obtained for a 
given degree of electron delocalization, i.e., for a given 
value of 7ij. For attractive interactions (C/ < 0) double 
occupations are favored and VI^[7jj] corresponds to the 
maximum of Ej(^iT"'U) ^or a given 7^ . 



B. Exact XC energy functional 

In order to determine £'c[7y] and II^[7ij] we look for 
the extremes of 






(1 -(*!*>) 

(2.12) 



with respect to |\I'). Lagrange multipliers e and \ij have 
been introduced to enforce the normalization of |^) and 
the conditions on the representability of 7^ . Derivation 



with respect to (^|, e and Ay yields the eigenvalue equa- 
tions 



(2.13) 



and the auxiliary conditions (^E'l^E') = 1 and 7y = 

(^1 Sct cJ^CjctI^)- The Lagrange multipliers Ay play 
the role of hopping integrals to be chosen in order that 
|\1/) yields the given 7ij. The pure-state representabil- 
ity of 7ij ensures that there is always a solution [[l9[ . 
In practice, however, one usually varies Ay in order to 
scan the domain of representability of 7y . For given 
Ay, the eigenstat e |^ o) corresponding to the lowest 



1^J = (*0lE, 
7., = (^lE.' 



eigenvalue of Eq. ( 2.13 ) yields the minimum W[7y] for 
cfcrCjal^o)- Any other I*) satisfying 
\^Cja\^) would have higher e and thus 
higher W . The subse t of 7 ^,- which are representable by 
a ground-state of Eq. (2.13) is the physically relevant one. 



since it necessarily includes the absolute minimum 7^ of 
El'jij]. Nevertheless, it should be noted that pure-state 
representable 7y may be considered that can only be 
represented by excited states or by linear combinations 
of eigenstates of Eq. ( 2.13| ). In the later case. Ay = Vij, 
and 1^0) is an eigenstate of the interaction term with low- 
est eigenvalue. Examples shall be di scussed in Sec. 
For the Hubbard model Eq. ( |2.13 ) reduces to 



[II 



E^. 



Cl^Ci 



I*) 



uY^i^^n 



*)==e|*>. (2.14) 



<ii> 



This eigenvalue problem can be solved numerically for 
clusters with different lattice structures and periodic 
boundary conditions. In this case we expand |^[7y]) 
in a complete set of basis states |$m) which have definite 
occupation numbers t'™ at all orbitals ia {niij\^m) — 
i/™|$™) with i/J^ = or 1). The values of u^^^ sat- 
isfy the usual conservation of the number of electrons 
Nfi = Nei + Nei and of the z component of the total spin 
5. = (TVeT - Nei)/2, where N^a = E^ C- For not too 
large cluster s, th e lowest energy |^o[7ij]) — the ground 
state of Eq. (2.14) — can be determined by sparse-matrix 



diagonalization procedures, for example, by using Lanc- 
zos iterative method Q. Iv&obij]) is calculated in the 
subspace of minimal Sz since this ensures that there are 
no a priori restrictions on the total spin S. In addi- 
tion, spin-projector operators may be used to study the 
dependence of £"0(712) on S. 

For a one-dimensional (ID) chain with NN hoppings 
ty — t, translational symmetry implies equal density- 
matrix elements 7y between NN's. T heref ore, one may 



set Ay = A for all NN ij, and then Eq. (2.14) has the same 



form as the ID Hubbard model for which Lieb and Wu's 
exact solution is available [g2|. In this case the lowest 
eigenvalue e is determined following the work by Shiba 
p3[. The coupled Bethe-Ansatz equations are solved as 



a function of A, band-filling n = Ne/Na, and for positive 
and negative U, by means of a simple iterative procedure. 



III. RESULTS AND DISCUSSION 

In this section we present and discuss exact results for 
the correlation energy functional Ec[jij] of the single- 
band Hubbard Hamiltonian with nearest neighbor hop- 
pings P|. Given the lattice structure, Na and Ne, the 
model is characterized by the dimensionless parameter 
U/t which measures the comp etitio n between kinetic and 
interaction energies [see Eq. ( |2.10|) ]. U > corresponds 
to the usual intra-atomic repulsive Coulomb interaction, 
while the attractive case {U < 0) simulates intra-atomic 
pairing of electrons. 



A. Repulsive interaction U > 

In Fig. I] the correlation energy Eq of the one- 
dimensional (ID) Hubbard model is shown for half-band 
filling (A^'e = Na) as a function of the density-matrix cl- 
ement or bond order 712 between NN's. 7y = 712 for all 
NN's i and j. Results are given for rings of finite length 
Na as well as for the infinite chain. Several general qual- 
itative features may be identified. First of all we observe 
that on bipartite lattices |^ £0(712) = £c(— 712), since 
the sign of the NN bond order can be changed with- 
out affecting the interaction energy 11^(712) by changing 
the phase of the local orbitals at one of the sublattices 
(ci^ —^ ~Cia for i £ A and Cj^ unchanged for j £ B, 
where A and B refer to the sublattices). Let us recall 
that the domain of definition of £0(712) is limited by 
the pure-state representability of 7y . The upper bound 
^12 ^^'^ the lower bound 7°^ for 712 (7°^ = —7°^ = 7?2 
on bipartite lattices) are the extreme values of the bond 
order between NN's on a given lattice and for given Na 
and Ne {■^ij — 712 for all NN ij). They represent the 



maximum degree of electron delocalization. 7j^2 ^'^d 
'-f-^^2 correspond to the extremes of the kinetic energy Ek 
[Ek = J2{ij) *y7y = izNa/2)t"fi2, where z is the co- 
ordination number] and thus to the ground state of the 
Hubbard model ior U — [7°^ for t > and 7°;^ for 
t < 0, see Eq. (2.10)]. For 712 = 712 the underlying elec- 
tronic state 1^0) is usually a single Slater determinant 
and therefore £0(712) = 0. In other words, the correla- 
tion energy vanishes as expected in the fully delocalized 
limit [g6|. As I712I decreases £0 decreases (£0 < 0) since 
correlations can reduce the Coulomb energy more and 
more efficiently as the electrons localize. £0 is minimum 
in the strongly correlated limit 712 = 7^ . For half-band 
filling this corresponds to a fully localized electronic state 
(71? = 0). Here, £0 cancels out the Hartree-Fock energy 
£hf and the Coulomb energy W vanishes (£g^ = — £hf) 
p?! . The ground-state values of 7f2 and Egg for a given 
iJ/t result from the competition between lowering £0 by 



decreasing 712 and lowering Ek by increasing it (t > 0). 
The divergence of dEc/dji2 for 712 = 712 is a neces- 
sary condition in order that 7f2 < 7i2 for arbitrary small 
U > 0. On the other side, for small 712, we observe 
that {Ec + ^hf) oc 7^2- This implies that for U/t > 1, 
7f2 oc t/U and Egg oc t'^ /U, a well known result in the 
Heisenberg limit of the Hubbard model (Ne — Na) JlO| . 

A more quantitative analysis of -^0(712) and in par- 
ticular the comparison of results for different Na is com- 
plicated by the size dependence of 7]'2 and E'hf- It is 
therefore useful to measure Ec in units of the Hartree- 
Fock energy and to bring the domains of representability 
to a common range by considering ec — Eq/Ey^f as a 
function of 512 = 712/712- Fig- H shows that £0(512) 
has approximately the same behavior for all considered 
Na- Finite size effects are small except for the very small 
sizes. The largest deviations from the common trend are 
found for Na = Ne — 4. Here we observe a discontin- 
uous drop of £c for 9i2 — 1 (512 < 1) which is due to 
the degeneracy of the single-particle spectrum. In fact 
in this case two of the four electrons occupy a doubly 
degenerate state in the uncorrelated limit and the mini- 
mum interaction energy W{^i2) does not correspond to 
a single-Slater-determinant state even for 712 = 7^2 [ p8[ . 
As Na increases £0(^12) approaches the infinite-length 
limit with alternations around the Na = c« curve. The 
strong similarity between £0(312) for small Na and for 
A^a = 00 is a remarkable result. It suggests that good 
approximations for £'0(712) in extended systems could 
be derived from finite cluster calculations. 

Fig. H shows the band-filling dependence of £'0(712) in 
a 10-site ID Hubbard ring. Results are given for N^ < 
Na, since for N^ > Na, £o(7*j,A^e) = £o(-7ii,2Afa - 
A'e) as a result of electron-hole symmetry ||2J]. Al- 
though £0(712) depends strongly on Ne, several qualita- 
tive properties are shared by all band fillings: (i) As in the 
half-filled band case, the domain of representability of 712 
is bound by the bond orders in the uncorrelated limits. 
In fact, 7°2" < 712 < 7i2^, where 7^*^ (712^) corresponds 
to the ground state of the U = tight-binding model for 
t > {t < 0). On bipartite lattices 7°^ = -7i2^ = 7i2- 
Notice that 712 increases monotonously with N^ as the 
single-particle band is filled up. This is an important 
contribution to the band-filling dependence of £0 (see 
Fig. |). (ii) In the delocahzed limit, £0(712) = for all 
the Ne for which W{'-f 12) derives from a single Slater de- 
terminant ||2q|. Moreover, the divergence of 9£o/97i2 
for 712 = 7f2 indicates that 7j'2 < 7^*2 for arbitrary 
small [/ > 0, as expected from perturbation theory, (iii) 
Starting from 712 — 712, £0(712) decreases with decreas- 
ing 712 reaching its lowest possible value E^ = — £hf 
for 712 = 7^"^ {Ne < Na). The same behavior is of 
course observed for 712 < 0. In particular, £0 = — £hf 
also for 712 = 7]^~- As shown in Fig. y, E^ decreases 
rapidly with increasing Ne, since £hf increases quadrat- 
ically with the electron density |2^]. (iv) On bipartite 
lattices 7^^"*" = — 7i2~ = 7i2i while on non-bipartite 
structures one generally has |7^^| / l7i^ L since the 



single-particle spectrum is different for positive and neg- 
ative energies. The decrease of £0 with decreasing I712I 
shows that the reduction of the Coulomb energy due to 
correlations is done at the expense of kinetic energy or 
electron delocalization, as already discussed for Ne — Na 
(Fig. p. (v) 7^^ > for all Ne < Na {j^ = for 
Ne = Na). 75I2 represents the largest NN bond order 
that can be constructed under the constraint of vanishing 
Coulomb repulsion energy. A lower bound for 7^^ is given 
by the bond order 7^^^^ in the fully-polarized ferromag- 
netic state (7]^ > '^12^'^)- This is obtained by occupying 
the lowest single-particle states with all electrons of the 
same spin {Ne < Na). Therefore, ^[2^'^ increases with A''e 
for Ne < Na/2 and then decreases for Na/2 < Ne < Na 
reaching 7^^*^ = at half-band filling {^[2^^ > for 
Ne < Na). In this way the non-monotonous dependence 
of 75^ on Ne can be explained (see Fig. g). (vi) The 
correlation energy is constant and equal to — £hf for 
7^~ < 712 < 7i2^- These values of 712 can never corre- 
spond to the ground-state energy of the Hubbard model, 
since in this range increasing 712 always lowers the kinetic 
energy {t > 0) without increasing the Coulomb repulsion 
(7S < 7f2 < 712)- For 7i^- < 712 < 7S+, 7i2 cannot 
be represented by a ground state of Eq. (2.14). In this 
range 712 can be derived from a linear combination of 
states having minimal Coulomb repulsion ||29|] . 

In order to compare the functional dependences of the 
correlation energy for different band fillings, it is useful 
to scale £0 in units of the Hartree-Fock energy and to 
bring the relevant domains 7J^ < 712 ^ 7i2 of differ- 
ent Ne to a common range. In Fig. ffl, £0 = £o/£hf 
is shown as a function of 512 = (712 — 712) /(7i2 ~ 7i?)- 
We observe that the results for £0(312) are remarkably 
similar for all band-fillings. The largest deviations from 
the common trend are found for Ne = 4. As already 
discussed for Na = Ae = 4, this anomalous behavior is 
related to the degeneracy of the single-particle spectrum 
and to the finite size of system. Fig. H shows that for 
the Hubbard model the largest part of the dependence 
of £0(712) on band filling comes from £hf, 7i2 and 7f^. 
Similar conclusions arc derived from the results for the 
infinite ID chain presented in Fig. H. For a given 312, 
£0(512) depends weakly on Ne/Na if the carrier density is 
low {Ne/Na < 0.4), and tends to increase as we approach 
half-band filling [see Fig. H(b)]. For high carrier densi- 
ties it become comparatively more difficult to minimize 
the Coulomb energy for a given degree of delocalization 
gi2. The effect is most pronounced for gu ~ 0.8-0.9, 
i.e., close to the uncorrelated limit. As we approach the 
strongly correlated limit {gu < 0.4) the dependence of 
£0 on Ne/Na is very weak even for Ne/Na ~ 1. One 
concludes that £0(512) is a useful basis for introducing 
practical approximations on more complex systems. 

The correlation energy £c is a universal functional of 
the complete single-particle density matrix 7^-. £o[7ij"] 
and Vl^[7ij] may depend on Na and Ne but are inde- 
pendent of tij and in particular of the lattice structure. 
The functional £0(712) considered in this paper depends 



by definition on the type of lattice, since the constraints 
imposed in the minimization only apply to NN bonds. 
In order to investigate this problem we have determined 
-6^0(712) for 2D and 3D finite clusters having Na < 12 
sites and periodic boundary conditions. In Fig. we 
compare these results with those of the ID 12-site pe- 
riodic ring. As shown in the inset figure, the quali- 
tative behavior is in all cases very similar. The main 
quantitative differences come from the domain of repre- 
sentability of 712, i.e., from the values of 7°^ and 7^2^ 
(7°^ < 7i2 < 7i^). Once scaled as a function of 712/7121 
Ec depends rather weakly on the lattice structure. No- 
tice that the Hartree-Fock energy E'hf = {U /A)Na is the 
same for all structures. However, for the BCC struc- 
ture we obtain VF(7i2) < i^HF, i.e., £'0(712) < 0; due to 
degeneracies in the single-particle spectrum of the con- 
sidered finite cluster [see inset Fig. ^(b)]. In order to 
correct for this finite size effect it is here more appropri- 
ate to consider ec = [^0(712) - £c(7i2)]/W"(7i2)- StiU, 
the differences in sq between BCC and FCC structures 
appear to be more important than between square and 
triangular 2D lattices. This is probably related to the de- 
generacies in the spectrum of the BCC cluster, as already 
observed for rings with Ne = 4m [Figs. and 0(a)]. 

The largest changes in ec for different lattice struc- 
tures are observed for intermediate degree of delocaliza- 
tion (512 ~ 0.7-0.9, see Fig. ||). Note that there is no 
monotonic trend as a function of the lattice dimension. 
For example, for gi2 = 0.7-0.9, ec first increases some- 
what as we go from ID to 2D lattices, but it then de- 
creases coming close to the ID curve for the 3D FCC 
lattice [ec(2D) > ec(FCC) > ec(lD) > ec(BCC) for 
0-7 < 512 < 0.9]. Finally, it is worth noting that in 
the strongly correlated hmit (512 < 0.3) the results for 
^0(512) are nearly the same for all considered lattice 
structures (see Fig. ^). This should be useful in order 
to develop simple general approximations to -Ec(7i2) in 
this limit. 



As in the repulsive case, 7^2^ < 712 < 712 since the 
domain of representability of 712 is independent of the 
form or type of the interaction. Moreover, -Ec(7i2) = 
_Ec("7i2) due to the electron-hole symmetry of bipartite 
lattices [|4|. Starting from 7°^ or 7^2^ (7°^ = —712' — 
7°2 on bipartite lattices), i?c(7i2) decreases with decreas- 
ing [712 1 reaching the minimum Eq" for 712 = 7^ and 
for 712 = 7i^~ (7i~+ = -7i^- = 7i^ in this case). For 
Ne even, M^(7i^) = iVeC//2, and for N^ odd, W{-^^) = 
(Ne — l)/7/2, which correspond to the maximum number 
of electron pairs that can be formed. For A^e even, the 
minimum E^ = U{Ne/2)[l- Ne/{2Na)] is achieved only 
for a complete electron localization (i.e, 7^ = 0). In con- 
trast, for odd Ne a finite-size effect is observed. In this 
case, one of the electrons remains unpaired even in the 
limit of strong electron correlations and the minimum of 



Ec is E; 



'^ = C/[(A^e-l)/2][l-(A^e + l)/(2A^a)]. Morcovcr, 
non- vanishing 7^^ are obtained as a result of the delocal- 
ization of the unpaired electron. 7^^ represents the max- 
imum bond order that can be obtained when (Ne — l)/2 
electron pairs are formed (7^ ^ for A'a ^ 00, A^e odd). 
Notice that in all cases the ground state 7f2 is found in 
the interval j^ < jf^ ^ 7i2- 

It is interesting to observe that i?c(7i2) can be appro- 
priately scaled in a similar way as for [/ > 0. In Fig. @(b), 
£0(^12) — Ec/\E^\ is shown as a function of the degree 
of delocalization gi2 = (712 - 7il)/(7i2 - 712)- £0(512) 
presents a pseudo-universal behavior in the sense that it 
depends weakly on A"a and Ne- The main deviations from 
the common trend are found for A'e = A^q = 4. As al- 
ready discussed for C/ > 0, this is a consequence of degen- 
eracies in the single-particle spectrum. In this case, the 
wave function corresponding to th e minimum in Levy's 
functional for 712 



7i°2 [Eq. (2.11 



cannot be described 
by a single Slater determinant and W{ji2 — > 712) < ^'hf- 



IV. CONCLUSION 



B. Attractive interaction [/ < 



The attractive Hubbard model describes itinerant elec- 
trons with local intra-atomic pairing ([/ < 0). The elec- 
tronic correlations are very different from those found in 
the repulsive case discussed so far. In particular Levy's 
interaction energy functional W{jij) now correspond to 
the maximum avera ge nu mber of double occupation for 
a given -fij [see Eq. (2.11)]. Therefore, it is very interest- 
ing to investigate the properties of the correlation energy 
functional Ecijij) also for U < and to contrast them 
with the results of the previous section. 

In Fig. the correlation energy -©0(712) of the attrac- 
tive Hubbard model is given at half-band filling for var- 
ious finite rings {Na < 12) and for the infinite ID chain 
(A'e = Na). The band-filling dependence of £^0(712) is 
shown in Fig. || for a finite 12-site ring (A'e < Na — 12). 



Density-matrix functional theory has been applied to 
lattice Hamiltonians taking the Hubbard model as a par- 
ticularly relevant example. In this framework the basic 
variable is the single-particle density matrix jij and the 
key unknown is the correlation energy functional Edjij]. 
The challenge is therefore to determine Ecljij] or to pro- 
vide with useful accurate approximations for it. In this 
paper we presented a systematic study of the functional 
dependence of £"0(712) on periodic lattices, where 712 
is the density-matrix element between nearest neighbors 
{lij — 712 for all NN ij). Based on finite-cluster ex- 
act diagonalizations and on the Bethe-Ansatz solution of 
the ID chain, we derived rigorous results for £0(712) of 
the Hubbard model as a function of the number of sites 
A"a, band filling Ne/Na and lattice structure. A basis 
for applications of density-matrix functional theory to 
many-body lattice models is thereby provided. The ob- 
served pseudo- universal behavior of eo(5i2) = £o/£hf 



as a function of 512 = (712 — 7i^)/(7i2 ~ 712) encourages 
transferring £0(512) from finite-size systems to infinite 
lattices or even to different lattice geometries. In fact, 
the exact i?c(7i2) of the Hubbard dimer has been re- 
cently used to infer a simple general Ansatz for -£^0(712) 
|pO| . With this approximation to i?c(7i2) the ground- 
state energies and charge-excitation gaps of ID and 2D 
lattices have been determined successfully in the whole 
range of U/t. Further investigations, for example, by 
considering magnetic impurity models or more complex 
multiband Hamiltonians, are certainly worthwhile. 
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FIG. 1. Correlation energy Ec of the Hubbard model on 
one- dimensional rings with Na sites and A^e ~ Na electrons 
as a function of the density-matrix element or bond order 
712 between nearest neighbors (NN). "ftj — 712 for all NN ij. 
U re fers t o the intra-atomic Coulomb repulsion [U > 0, see 

Eq. (ETol)] pa. 
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FIG. 2. Correlation energy Ec in units of the 

Hartree-Fock energy _Bhf [^ for the Hubbard model on 
one- dimensional rings. Results are given as a function of 
712/712, where 712 stands for the NN ground-state bond order 
in the uncorrelated limit {U — 0). jij — 712 for all NN ij. Na 
refers to the number of sites, and Ne — Na to the number of 
electrons. £0(712) = -Eo(-7i2), see Fig. |l|. 
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FIG. 3. Band-filling dependence of the correlation 

energy £0(712) of the one-dimensional Hubbard model 
for Na = 10 sites. ^ij = 712 for all NN ij. 

Ne refers to the number of electrons and U to the 
intra-atomic Coulomb repulsion. On bipartite lattices 
Ec{'yi2,Ne) = Ec{-'yi2,N,) = £0(712, 2iVa - N,) pi. 
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FIG. 4. Correlation energy Ec in units of the 

Hartree-Fock energy _Ehf |^ for the one-dimensional 
Hubbard model on a ID 10-site ring. Results are 
given as a function of the degree of delocalization 
ffi2 = (712 -712)/ (712 —712), where 7^2 refers to the NN bond 
order in the uncorrelated ground-state {U — 0) and 712 to the 
NN bond order in the strongly correlated limit {U/t —f 00). 
As in Fig. pi difi'erent band-fillings Ne/Na are considered: (a) 
iVe < 6 and (b) iVe > 6. 
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FIG. 6. Correlation energy £^0(712) of the Hubbard model 
on different lattice structures. Finite clusters with periodic 
boundary conditions are considered at half-band filling: (a) 
one- dimensional (ID) ring (Na — 12), 2D square and trian- 
gular lattices (3x4 clusters), and (b) 3D fee and bcc lattices 
(two-tetrahedron cluster with Na = 8) Pq|. U refers to the 
intra-atomic Coulomb repulsion ([/ > 0). Notice the effect of 
scaling 712 with the uncorrelated 7^2 by comparing main and 
inset figures. 



FIG. 5. Correlation energy Ec of the Hubbard model 
on the infinite one-dimensional chain. Results are given for 
£c ~ Ec/EuF (-Ehf = Hartree-Fock energy) as a function of 

(a) gi2 = (712 - 7SL/(7?2 - 7S), and (b) band filling N^/Na 
(see caption of Fig. m . 
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FIG. 7. Correlation energy Ec of the attractive Hubbard 
model ([/ < 0) on one-dimensional rings with Na sites and 
A^e ~ Na electrons. In (a) Ec is shown as a function of the 
density-matrix element 712 between nearest neighbors (NN), 
-yij = 712 for all NN ij. In (b) Ec/\Eq'\ is given as a function 
of the degree of delocalization 1712 = (712 — 712)/ (712 ~ 712)- 
See caption of Fig. W. 
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!. Band-filling dependence of the correlation energy 
of the one-dimensional attractive Hubbard model 
{U < 0). The number of sites is Na = 12 and the num- 
ber of electrons A^e are indicated. In (a) Ec is shown as a 
function of 712 and in (b) Ec/lEc'l is given as a function 
of the degree of delocalization gi2 = (712 — 7i2)/(7ii,~ 7S) 
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